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Abstract 

The fractional Poisson process has recently attracted experts from several fields of study. Its 
natural generalization of the ordinary Poisson process made the model more appealing for real-world 
applications. In this paper, we generalized the standard and fractional Poisson processes through 
the waiting time distribution, and showed their relations to an integral operator with a generalized 
Mittag-Leffler function in the kernel. The waiting times of the proposed renewal processes have 
the generalized Mittag-Leffler and stretched-squashed Mittag-Leffler distributions. Note that the 
generalizations naturally provide greater flexibility in modeling real-life renewal processes. Algorithms 
to simulate sample paths and to estimate the model parameters are derived. Note also that these 
procedures are necessary to make these models more usable in practice. State probabilities and other 
qualitative or quantitative features of the models are also discussed. 

Keywords: Fractional Poisson process, generalized Mittag-Leffler distribution, renewal processes, 
Prabhakar operator. 



1 Introduction 

The fractional Poisson process HT U13II gained popularity in many areas of research as it naturally 
generalizes the standard or classical Poisson process. Recall that the inter-event time density function of 
the fractional Poisson process N v (f), r > 0, v e (0, 1], was originally derived in Repin and Saichev ||14ll 
(known to date) and has the following integral form: 



where 



e" x (/> v (At/x)dx, ve(0,l], t >0,A>0, (1.1) 



sin(v7i) 

<M£>= r;>± ' , — Ci.2) 

rc[s +s +2cos(v7i)] 



The preceding density function suggests that the tail distribution of the waiting time is of the form 



Pr(r>t) = £ v (-At v ), 



(1.3) 



where 

EpW = 2] zeC,/3eC,8l(/3)>0 ) (1.4) 

is the Mittag-Leffler function. Note that the Mittag-Leffler density has been widely used to describe 
distributions appearing in anomalous diffusion, finance and economics, transport of charge carriers in 
semiconductors, and light propagation through random media (see, e.g., II15II16II "). In view of equations 
( |1.3D and Q1.4D , the interarrival time density for the fractional Poisson process directly follows as 

f v {t) = Xt v - l E ViV {-kt v ), t>0, (1.5) 



where 



E. 



3 , r (z) = V / zeC, 0,reC, 8l(j8)>0 (1.6) 

is the two-parameter Mittag-Leffler function. The qth fractional moment ||17ll of the random interarrival 
time is 

7ir(l + q) 

E [ r T = 3«,rf TW^r 7 W1 V 0<q<v. (1.7) 

A 9 r(q/vjsin(7rq/vjr(.l — qj 

In addition, the above information automatically gives the probability density function 

f.vm—1 

/ - (0 = Am o^iji^" 1)( - AtV) ' (L8) 

of the m-th arrival time because its Laplace transform, 



^(Odt=^^, (1.9) 



where At v ) is the fcth derivative of E vv (z) evaluated at z = —At 1 '. As v — » 1, the above 

distribution converges to the classical Erlang distribution. 

In another approach to the study the fractional Poisson process, Laskin [lj used the fractional 
Kolmogorov-Feller-type differential equation system to characterize the one-dimensional state probabil- 
ity distributions as (see Laskin [1, formula (25)] and Beghin and Orsingher [5, formula (2.5)]) 

(At v ) fc ~ (r + fc)! (-At r ) r 



r=0 

One can also show [1 ] that the moment generating function (MGF) of the fractional Poisson process is 



pl(t) = MN v (t) = k}=^— -> ' - ,. , n , fc>0, t>0. (1.10) 

fc! ^—i r! r(v(r + fc)+lj 



^\Xt v (e- s -l)Y 
M v (s, t) = E v [Ke~ s - l)t v = ^ L r ; ^./ J , (LID 

n r(vr + l) 

which permits calculation (see Table [T]) of the moments. A summary of the characteristics of the classical 
and fractional Poisson processes is shown in Table [T] below. 

In this paper, we generalize the standard and fractional Poisson processes through their waiting 
time distributions. In particular, we propose two renewal processes that have waiting times that are 
generalized Mittag-Leffler and stretched-squashed Mittag-Leffler distributed. These generalizations 
naturally provide more flexibility in capturing real-world renewal processes. Algorithms to simulate 
sample paths and estimate the model parameters are derived and tested. State probabilities and other 
qualitative or quantitative features of the models are also discussed. 

The rest of the paper is organized as follows. In Section |2j a renewal process with generalized 
Mittag-Leffler distributed waiting times is presented. Procedures to generate sample paths and to 
estimate parameters are also derived. In Section [3] another generalization based on stretching and 
squashing the Mittag-Leffler distributed inter-event times is developed. Methods to simulate sample 
trajectories and to estimate parameters are also showcased. More discussions are provided in Section [4] 
Finally, computational test results are shown in the appendix. 
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Poisson process (v = 1) Fractional Poisson Process (v < 1) 



Pr(T v > t) 
/ v (0 

Mean 
Variance 
fcth moment 



e~ Xt 
Ae" At 



At 
At 



E v (-Ar v ) 
A,t v_1 £ VjV (-At v ) 

{Xff yoo (r+fe)! (-At v ) r 

At v /r(v + i) 



A t 

r(v 



(-l)^exp[A(e- 5 -l)t] 



5=0 



•+i) + ^ [ vr(2v) r 2 (r+i) ] 



Table 1 : Properties of fractional Poisson process compared with those of the standard Poisson process. 

2 Generalization I 

We consider the generalized Mittag-Leffler distribution (see e.g. Pillai [18]) built from the generalized 
Mittag-Leffler function ||19H20|| . Let T v ' 5 be a generalized Mittag-Leffler distributed random variable. 
Then the probability density function is 



f v -°(t) = A 5 t 5v - 1 E v 5 j5v (-At v ), t > 0, A > 0, v e (0, 1], 8 . 



where 



4rW = S 



(a 



^r!r(/3r + r ) 



(2.1) 



(2.2) 



is the generalized Mittag-Leffler function (see Figure]!]). The Pochhammer symbol (£) r can be written 
also as (£) r = + 1) . . . (£ + r - 1), £ ^ 0. When 5v < 1 the function ( |2.lD has an asymptote at t = 0, 
while in the particular case 5v = 1 



/ v ' 5 (t)| £=0 = A 1 ^ 1 v (-At v ) 
The Laplace transform of ( |2.lD reads 

L{A 5 t 5, '- 1 £ t % r (-At 1 ')}(s) = 



: A 



1/v 



(s v + A) 5 



(2.3) 



(2.4) 



(see Mathai and Haubold [21], formula (2.3.24), page 95). Below are the plots of the generalized 
Mittag-Leffler densities. 

Now, consider m i.i.d. random waiting times ^, i = l,...,m of a renewal point process, here 

denoted as iV v ' 5 (t), t > 0, which are distributed as in ( |2.1| ). Furthermore, denote V n ' 5 = <? x H Y3" m 

as the waiting time of the mth renewal event and T^' 5 = 0. Then 



Ee~ 



A 5 



(s v + A) 5m ' 

By formula (2.3.24) of Mathai and Haubold [EU we have 

Pr{:r< 5 e dt}/dt = X 5m t v5m - l E 5 ™ 5 J-Xt v ), t > 0, A > 0. 



(2.5) 



(2.6) 



3 




Figure 1: The generalized Mittag-Leffler density plots (see ( |2.1| )) for parameter values (5, A) = (0.5, 1) 
(top) and (5, X) = (2, 1) (bottom) where v goes from 0.2 to 1 with step size of 0.2. 
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It is rather immediate now to obtain the state probabilities pi' (t) = Pr{JV v ' 5 (t) = k}, k > because 



e" st Pr{iV v - 5 (t) = fc}dt 



(2.7) 



o 

f CO 



e -st (p r{r v.5 <t} _p r{r v,5 <t ^ dt 





f CO 



Pr{T, v ' 5 edy}- 



PrlT^edyl^dt 



Pr{r- 6 edy} 



e" st dt- 





00 



e" st dt 



e - S yp r{r v,5 edy} _ 



A 



5fc 



fc> 0. 



Inverting the preceding Laplace transform we readily arrive at 

p v,5 (0 = ^k^k^^^ _ A5 (, + l) t v 5 (, + l) £ ^l) +i)+i( _ At1); k > _ (2 8) 

When (0) r = 0, r e N, (0) = 1, and k = 0, equation ( |2.8| ) becomes 

P ;'- 5 (t) = i-A 5 r 5 £ i 5 r5+1 (-At 1 '). 

Clearly, 

1, k = 0, 



(2.9) 



p, v ' 5 (0) = Pr{AT- fl (0) = k} 



r,5f 



o, fc>i. 



(2.10) 



Theorem 2.1. The state probabilities pi' (t) = Pr{JV v ' 5 (t) = fc}, fc > 1, satisfy the convolution-type 
Volterra equation of the first kind 



/o 



v,5 



(2.11) 



Proof. We start by rewriting fl2.7| ) by means of the following well-known relation (see e.g. Haubold 
et al. [22], formula (11.7), page 17): 



f (x - tf-^ [a(x - t) a } t v - 1 B^ v (at 8 )dt = x^-^J^ax*), 
Jo 

where a, /3, 7, a, v, cr e C, and 31(a) > 0, 3t(/3) > 0, 31(f) > 0, 3t(v) > 0, 2H(cr) > 0. Then 
pl'\t) = X 5 h v5k E 5 v % k+x {-Xt v ) - l S ™t vB <* + »E^ +l ^-Xn 



(2.12) 



(2.13) 



A 



(t -5) v5(fc -XvS-iHi [_ACt -*) v ]* vff -Xv5C-^ v )ds 



r f 



. A 5(Jc+l) 



(t -s) v5fc E v % +1 [-A(t -s7]s ve - 1 £ v % 5 (-^ v )ds 



w=t-s ^5 



(t - wr 5 - 1 ^,,, [-A(t - w) 1 '] A^V^^K- ( _ Aw v )dw 



-A 5 f {t-wy 5 - l E 5 vv5 [-Kt-wy]X 5k w v5 *Ell 5k+l {-Xw v )&w 
Jo 



(t-w) v5 - l E v 5 v5 [-A(t-w) 1 '] 



[a^^V 



5(it-l)p5(/c-l) 



x 8k w v5k E 5k 



J v,vS(fc-l)+l 



v,v5fe+l 



(-Aw 1 )] 



dw 



r f 



(t - wf-^ [(-A(t - wDlp^CuOdw. 



□ 



Remark 2.1. Result ( |2.11| ) can aZso be conveniently expressed by means of the Prabhakar operator [20], 
defined as 

{K,y.^,aA) M = (x - y) M_1 ^ OO - yY) 4(y)dy, x > a, p,[i,y e C, 9t(p), R(/i) > 0, 

J a 

(2.14) 

which is a generalization of the Riemann-Liouville fractional integral of <p(x). Therefore, we obtain 

pl'\t) = X 5 {E 5 vv5 _ x . Q+ pl^)M. (2.15) 
Remark 2.2. When 5 = 1 equation fl2.ll! ) clearly reduces to 

(2.16) 



p£(t) = A (t-w7- i £ v , v (-A(t-wr)p'_ 1 (w)dw. 



We now check that the state probabilities pj. (t) o/ the fractional Poisson process JV v (t), t > (see e.g. 
Beghin and Orsingher /|5l7 ) satisfy the above integral equation. By recalling that 



p2C0 = (At v )X+ fc+ i(- AtV )> fc > o, t > 0, 



we can wnte 



Pi'Ct) = A (t - MO v_1 Bv,v(-A(t - dw 
Jo 

= A (t - wf-^^C-ACt - w7) (Aw) fc - 1 £ v , v(fc _ 1)+1 (-Aw 1 ) dw 
Jo 

Notice also that for v = 5 = 1 (classical case), equation fl2.11| ) reduces to 

r f 



(2.17) 



(2.18) 



-A(t-w) 



p fe -i(w)dw, 



(2.19) 



where Pk(0, k> 1, t > 0, are the state probabilities of a homogeneous Poisson process JV(t), t > 0. 
Equation fl2.19| ) is easiZy solvable and the solution reads 



p k -i(t) = X- 1 — Pk (t) + p k (t), k>l. 
at 



(2.20) 



Finally we note that the above equation is clearly the difference- differential equation governing the state 
probabilities of a homogeneous Poisson process. 
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Theorem 2.2. The state probabilities pi ' (t) = Pr{JV r,5 (t) = k], k > 0, t > 0, satisfy the equations 



d 1 



5+0 



dt 



'5+e 



( E v"l-A;o + Pl' 5 ) (0 = ^P&CO + 5 W (t- v5 E-f_ v5 (-At v ) - A 5 ) , (2.21) 



/or any 9 e C, 31(0) > 0, where 5 kfi is the Kronecker's delta and where the operator £ vS +e ^ tae 
Riemann-Liouville fractional derivative of order v5 + 9. 

Proof. We start by considering k > 1. Applying the the operator E~ 5 e _ A . 0+ to equation ( |2.15| ), we obtain 

tei-AiW-Pfc 6 ) CO = ^ (e;1_ A;0+ « v5 ,-. ; o + P fc V l 5 i)) CO (2-22) 

~ fclw* 5 ) co = a^p^co, 

where J t V o+ e is the Riemann-Liouville fractional integral operator. By recalling that the Riemann- 
Liouville fractional derivative is the left inverse operator of the Riemann-Liouville fractional integral 
(see e.g. Diethelm [23], Theorem 2.14, page 30), we readily arrive at the claimed result. For k = 0, it is 
sufficient to show that 



d v 5+ e _ 

dr 5+e 



•' Pi' 5 ) CO 



-AjO+i'O 



(2.23) 



d r5+0 

dr 5+9 

AvS+e 

dt v5+e 
r t 

-A 5 



(t - yf-%1 [-A(t - yf ] p^(y)dy 



(t-y) 6 "^ [-A(t-y) 1 '] dy 



(t-y) 9 " 1 ^ [-A(t-y) v ]y v5 E v 5 v5+1 (-Ay v )dy 

d v5+e /- n 
= nzs I t e E _ * (-At v )-A 5 t v6+e £° ^ fli1 (-At v )) 



5+0 



d v 



5+0 



dt 1 



v5+e 



r(v5 + + l) 



r v5 E~? J-Xt r )-X°. 

v.l— vo v J 



□ 



For more information on the inverse operator appearing in ( |2.2lD , the reader can consult Saigo et al. 
191 , Section 6. 



Remark 2.3. When 5 = 1 and k>l, equation Q2.21| ) can be written as 

d r+f 



pl_ 1 {t) = X~ 1 



dt r+0 
d r + 



(t - wy-'E-X-Xit - w? M(t)dw 



(2.24) 



dt 1 ' 



r+0 



(t-w) 6 " 1 



1 1 

+ A(£ - w) v 



r(0) 



r(v + 9) 



p?'(w)dw 



= A" 



L d v d e 1 

d^dt^rcej 



Ct-w)°-^(w)dw + -— 



dt v+0 r(v + e) 



(t-wr +e -V(w)dw, 



while, for k = 0, and considering that 

r v E~] (-At v )-A=r v 

v.l— v v J 



1 A t 

+ At 1 ' - A 



r(i-v) 



r(i-v)' 



(2.25) 
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we have 

r' 



jv+9 r t 



T(l - v) dt v+e 

d v+9 



(t - w) e - 1 £; 1 9 (-A(t - w7) P ;(t) dw 



(2.26) 



dt v+e 
d v d 6 



(t-w) 9 " 1 



1 1 

+ A(t-w) v 



r(0) 



dt v dt e r(0) 



r(v + e) 

d v+9 



Po'(w)dw 



(t-wrxwd W +A- ^ 



dt v+e r(v + 0) 



(t-w) r+9 - 1 pl'(w)dw. 



Hence, we retrieve the fractional difference- differential equations governing the state probabilities of a 
fractional Poisson process £T]/: 



—PIM = -Xpl{t) + XpU{t) + 5 k , 0f ^- y k > 0. 



(2.27) 



From equation Q2.21| ) we can easily arrive at the following partial differential equation for the 
probability generating function % 5 (u, t) = Y^kLo^Pk' 5 ^- 



3v5+e 



(e;1_ A; o + ^, 5 (" 5 O) (0 = A 5 n «^ s (u, t) + r v5 E-f_ v5 (-At v ) - A 5 . (2.28) 

From the above equation and by recalling the formula J^ Vj5 (u, t)\ u=1 = EJV 1,5 (t), it is now immediate 
to derive the differential equation involving the mean value as 



d r 



r5+e 



dr 



^ (E-^_ A;O+ EiV^ 5 (0) (0 = A 5 (l +EI\T' 5 (t)) . 



(2.29) 



Observe that equations ( |2.28| ) and ( |2.29| ) reduce to the corresponding equations in the pure fractional 
case when 5 = 1 (see Laskin [ 1 , formula (22)] for the differential equation involving the probability 
generating function). For the fractional Poisson process N v [t), t > (5 = 1), equation ( |2.29| ) becomes 



d v+e 

dt r + ^ E v,e,-A;0+ 
d r+0 r t 



(e^.^o+EATG)) (t) = A + AEAT(0 

(t - yf- l E-\ [-A(t - y) v ] EAT(t) dy = A + AEAT(t) 



(2.30) 



r+6 



dt v 
d v d 



9 



r f 



dt v dt e r(0) 

d r+9 



+ A 



dt v+e r( - v + ) 



(t-y) e - 1 EAT v (y)dy 
r f 

(t - y) v+e_1 EiV r (y)dy = A + AEAT(t) 



d v 



dt 



-EAP(0 = A, 



with EiV v (0) = and considering that the second step is justified by the semigroup property of the 
Riemann-Liouville fractional derivative (see Diethelm [23, Theorem 2.2, page 14]). The solution to 
H2.30D is well-known and reads [1, formula (26)] 



EAP(t): 



At 1 



v e (0, 1] t > 0. 



nv + i)' 

The following theorem derives the mean value of the process N v,5 (t), t > 0. 



(2.31) 
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Theorem 2.3. Let v e (0, l],5eC,9e C, 31(0) > 0. The solution to 



dt* 



(E; 5 9 _ A . 0+ EN V " 5 (-)) (t) = A 5 (l + EAP- 5 (t)) , 



[ gSS ( e ;1-a;o+ e N v,5 (0) (£)] t ^ Q = 0, Vfc = 0, . . . , n - 1, n - 1 < 8t(v5 + 9) < n, 



reads 



EN 



- 5 (t) = £A 5 ( r+ «t 1 ' 5 ^ +1 ^^ 5 + ( 1 r ) +1)+1 (-At 1 '). 



Proof. We start by taking the Laplace transform of ( |2.32| ), obtaining 



v5+9 



dt 



■5+8 



(E; 5 9| _ A;O+ EiV- 5 (0)(t)dt = ^ + A 5 



<^> s 



v5+B 



v8+9 



r f 





f oo 



e- st EN r -°(t)dt 
A 5 



(2.32) 



(2.33) 



(2.34) 



(t-y) 9 " 1 ^ [-A(t-y) r ]EAP' 5 (y)dydt = — + A 5 \ e~ 5t EN v ' 5 {t)dt 



EN v ' 5 (y)dy 



o 

f CO 



'5+e 



y 

f CO 



tft-yf-^-j [_A(t-yr] dt = — + A S 



e- 5t EJV r - 5 (t)dt 



EJV v > 5 (y)dy 



o 

/- 00 



e -*C«+^ z e-i£-«(_Aj[ v )dj[ = — + A a 



e - st MN v - 5 (t)dt 



e - 5J, EiV v ' 5 (y)dy 



e - S z z e-i E -5 Az= _ + X 5 
v,e s 



e - st EN v '"{t)dt 



<=> s v5+e L {EAP'' 5 (t)} (s)s" e (l + As" 1 ) 5 = s^A 5 + A 5 L {EiV r - 5 (t)} (s) 

A 5 

<=> L{EJV v ' 5 (t)l(s)= — - — . 

Notice that the first step in Q2.34| ) is justified by the formula for the Laplace transform of the Riemann- 
Liouville fractional derivative and by applying the initial conditions. Before inverting the Laplace 
transform, it can be shown that 



L{EAT' 5 (t)}(s): 



A 



(2.35) 



= s(s v + A) 5 ^ 

-! 00 

c £—1 i 



A 



(s v + A) a 



A 5(r+1) 



Thus, the mean value is now easily found by inverting ( |2.35| ) term by term: 

00 ft 

EAP< 5 (t) = J] A^ + V 5(r+1) XvS?r + i)(-^ V )^ 

r=0 Jo 

oo 

= y A 5(r + l)v5(r+l) £ 5(r+l) ( _A t V). 
/ i v,vo(r+l)+l v - y 



(2.36) 



r=0 



□ 
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Remark 2.4. For 5 = 1, the mean value ( |2.33| ) reduces to that of the pure fractional case ( |2.31| ). Indeed, 



r=0 

and passing now to the Laplace transform we get 

00 

Z.{EJV V (0}GO = ^A r+1 s-'' (r+1) (1 + As- r )" (r+1) = X/s' 



r=0 



which immediately leads to ( |2.31[ ). When 5=1, ( |2.32| l reduces to 
jv+e-Jc-i 



dt 



M-8-k- 



: 



(2.37) 



(2.38) 



(2.39) 



JV+0-fc-l 

jjj-v+e-fc-i 



r f 



o 

r f 



(t-y) 9 " 1 



1 l 

+ A(t-y)' 



rce) 



r(i + 0) 



EAT(y)dy 



— EAP'(t) + AEAT(t) 



0, 



/or each fc = 0,...,n-l J n-l< 3l(v + 6) < n. By recalling E«/K(t) = £™ rpj^(t) and equation ( |2.10D 
/or 5 = 1 we obtain 



d r-i 

-EAP'(t) 



0, 



(2.40) 



where we considered only k = and therefore only one initial condition is used. 



2.1 Path simulation and parameter estimation 

It is straightforward to generate a sample trajectory of generalization I by noting that the generalized 
Mittag-leffler random variable T v,s (see, e.g., Pillai 111811 ) is a mixture of gamma densities, i.e., 



T v ' = U 1/V V V , 
where U is gamma distributed with density function 

X 5 

f u (u)=—u 5 - 1 e- M , u>0, 
HP) 



(2.41) 



(2.42) 



and V v is strictly positive-stable distributed with exp(— s r ) as the Laplace transform of the corresponding 
density function. Note that the qth fractional moment of the inter-event time can be easily shown as 



511 



e [r- 5 ] 



nT{q/v + 5) 



A9/ v r(q/v)sin(7iq/v)r(l-q)' 



< q < v. 



(2.43) 



Typically, generating T v,s and adding one (corresponding to a single jump or event) each time gives a 
sample trajectory. 

Given m jumps (corresponding to m renewal times), we propose method-of-moments estimators 
for the parameters v, 5, and X to make the preceding generalization usable in practice. Getting the 
logarithm of T v > 5 we have 

T'=-U' + V' (2.44) 
v 
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where T = ln(T v ' 5 ), U' = ln(LO, and V' = ln(V v ). Following Cahoy et al. [6] we get the estimating 
equations: 

(\ \ t/>(5) — ln(A) 
(i P = E(T')=» } — 1 +— — , (2.45) 



v 

^2) (5) _ 2 (v 3 -l)C(3) 
ju 3 =E(r / -Mr) = S . ( 2 -47) 

y] » 0.57721 is the Euler's constant, and £(3) is the Pdemann Zeta function evaluated at 3. Using the 
equations of the variance and the third central moment above, we can solve for the estimates 5 and v 
using pL 3 and &j, . Plugging v and 5 into the mean equation above, we obtain the estimate of A as 

A = exp (- [v (At' - fiW ~ 1)) - tK$)] ) • (2-48) 

Furthermore, we tested the above procedure using the following estimate of the digamma function: 

i/»(t) = log(T) - 1/(2t) - 1/Q2t 2 ) + 1/(120t 4 ) - 1/(252t 6 ) + 0(1/t 8 ). (2.49) 

We then calculated the bias and the root-mean-square-error (RMSE) based on the 1000 generated data 
samples for different parameter values. Table [2] in the appendix generally indicated positive results for 
the proposed method. 

3 Generalization II 

Recall that a random variable X is Mittag-Leffler-distributed with parameters A > and v e (0, 1] if it 
has probability density function 

f x {x) = Xx v - 1 E ViV {-Xx v ), xeK+, (3.1) 

where 

Ea,pW = Y. rf * v xeK, (3.2) 

is the Mittag-Leffler function. Note that Vr{X > x] = E v l {-\x v ). 

Let Y = 1/X. Then the random variable Y has the inverse Mittag-Leffler distribution, that is, 

Pr{Y <y}= Pr{l/X < y} = Pr{X > 1/y} = E Vil {-Xy~ v ). (3.3) 

Hence, the corresponding probability density function is 

— Pr{7 < y} = -E v ^-Xy-^vy-^ 1 (3.4) 
dy v 

= Xy- v - 1 E ViV (-ly- v ), jeR + . 

When v = 1, formula ( |3.4P is the probability density function of an inverse exponential random variable, 
that is, 

f yeK + . (3.5) 

y 
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We now give a single definition for both the Mittag-Leffler and the inverse Mittag-Leffler distribu- 
tions. Note that the probability density function in equation ( |3.3| ) can be written as 



/y(y) = Ay r_1 V(-Ay r ), y e R + , v e (0, l], (3.6) 
where 7 = ±v. By freeing the parameter 7 in formula ( |3.6| ), we arrive at the probability density 

Irl, 



/ B (C) = m^XvC-Ar), ? e K+, v e (0, 1], 7 e R\{0} 
v 

(see Figure^. Observe that J™ f s (E,)dE, = 1 as 

I I f CO 1 1 /" CO 

Irl, 1 kv_i„ , ^c?=^ v/r ) Irl. 



v 

/* CO 



Irl 



* v - v/ % v (-A« v )* v/ r- 1 dai 



(3.7) 



(3.8) 



Az , '- 1 £ vv (-Az r )dz= 1. 



Note that in the second-to-last line of ( |3.8P , sgn(y) is used to stabilize the domain of the integral. 
The Laplace transform Ee~ s ^ can be shown as 



Ee" sC = A 



Irl 



(3.9) 



■2^1 



-As" v 



J] (-As" 1 ') 



(1,1), (r,v) 

(v,v) 
r T(rv + 7) 



- r=0 r(rv + v)' 

where we use formula (2.2.22) of Mathai and Haubold [12111 . When 7 = v e (0, 1] we obtain 

A 



Ee" 



s v + A' 



(3.10) 



as in Beghin and Orsingher O, formula (4.15). 

From the above discussion it is clear that S = X vlr , v e (0, 1], 7 e R\{0}, where X is the Mittag- 
Leffler distribution with probability density function ( |3.1[ ), therefore the waiting times S of a newly 
constructed renewal process are simple time-stretching or time-squashing of the original Mittag-Leffler 
distributed waiting times X. Below are density plots of S = X vlr where X has the generalized Mittag- 
Leffler distribution given in ( |2.1| ). 

Now, consider m i.i.d. random inter-event imes E 1 ,...,E m of a counting process jY{t), t > 0, which 
are distributed according to (|3.7D- If W m = E 1 H h 3 m is the waiting time till the mth event, we have 



Ee- sW '» =s- mr I A- 



2^ 



-As" v 



(1,1), Cr.v) 

(v,v) 



To determine the state probabilities q£(t) = T?x{jY(t) = k},k> 0, we write 



e- st qUt)dt 



e~ st (Pr{W k < r} - MW k+1 < t}) dt 



(3.11) 



(3.12) 



-As" 



(1,1), Cr.v) 

(v,v) 



-(fc+iJr-i A 



Irl 



. k+1 




) 2^ +1 


-As" 1 ' 



(i,i), Cr.v) 

(v,v) 
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Figure 2: The stretched-squashed Mittag-Leffler density plots (see ( |3.7| )) for parameter val- 
ues (v,8,X,r) = ((0.1, 0.5,1), 2, 1,-0.5) (top), (v,5,A,y) = ((0.1,0.5, 1), 2, 1, 1) (middle), and 
(v, 5, A, y) = ((0.1, 0.5, 1), 2, 1, 5) (bottom). 
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3.1 Path generation and parameter estimation 



Simulating a sample path of generalization II directly follows from generalization I. Note that the S's 
can be generated using the algorithm of Cahoy et al. [6! • It is also straightforward to show that the qth 
fractional moment of the random inter-event time is 



7ir(qv Ij + 1) 



A« v ^r(q/v)sin(7iq/v)r(l -q)' 



< q < v. 



(3.13) 



Given m renewal times, we propose a formal procedure to estimate the parameters v, j, and A of 
generalization II. Let S' = ln(S) and X' = ln(X). Following Cahoy et al. [6], we can deduce that 



Me 



:E(S') 



and 



/i 3 = E(S / -/ia')' 



-ln(A) 

v 

i : 
3v 2 < 

-2CO) 



(3.14) 



(3.15) 



(3.16) 



Using the estimating equations above, we can eliminate y and solve for v by getting the 2/3 root of the 
third central moment and dividing it by the variance. Thus, we obtain 



^ 3 [(2C(3)) 2/3 + f ]' 
where c = {pfj 3 ) / . Substituting v to the variance equation fl3.15| ), we get 



(3.17) 




1 

3^2 



Finally, plugging v and f into the mean equation ( [3.14D above, we have 

X = exp [- (As'f + 17*)] . 



(3.18) 



(3.19) 



We also tested the above explicit forms of the estimators by calculating the bias and the root-mean- 
square-error (RMSE) based on the 1000 generated data samples for different parameter and total jump 
size values. Overall, Table|3]in the appendix showed favorable results for the proposed procedure. 



4 Concluding remarks 

We proposed two generalizations of the standard and the fractional Poisson processes through their 
renewal time distributions which naturally provided greater flexibility in modeling real-life renewal 
processes. Statistical properties such as the state probabilities and process moments were derived. 
Algorithms to simulate trajectories and to estimate model parameters were also developed. Generally, 
tests provided additional merits to the proposed procedures. 

Although some work have already been done, there are still a few things that need to be pursued. 
For instance, the complete analysis of the counting process related to the renewal process that has 
stretched-squashed generalized Mittag-Leffler distributed waiting times would be a worthy pursuit. 
Also, the development of estimators using likelihood approaches would be of interest as well. 
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5 Appendix 









Bias 






RMSE 




V 


Est 


m= 100 


1000 


10000 


m = 100 


1000 


10000 




V 


0.018 


0.001 


0.000 


0.184 


0.049 


0.013 


0.5 


5 


0.082 


0.011 


0.002 


0.253 


0.074 


0.023 




X 


0.207 


0.027 


0.003 


0.525 


0.151 


0.048 




V 


0.026 


0.003 


0.000 


0.210 


0.068 


0.016 


0.6 


5 


0.065 


0.010 


0.000 


0.205 


0.074 


0.023 




X 


0.174 


0.026 


0.001 


0.442 


0.151 


0.047 




V 


0.024 


0.001 


0.000 


0.254 


0.056 


0.018 


0.7 


8 


0.069 


0.009 


0.000 


0.207 


0.067 


0.023 




X 


0.176 


0.022 


0.001 


0.434 


0.135 


0.046 




V 


0.005 


0.003 


0.000 


0.264 


0.077 


0.020 


0.8 


8 


0.080 


0.008 


0.001 


0.199 


0.072 


0.023 




X 


0.202 


0.022 


0.004 


0.429 


0.147 


0.047 




V 


0.022 


0.001 


0.000 


0.349 


0.079 


0.021 


0.95 


8 


0.070 


0.009 


0.002 


0.187 


0.067 


0.022 




X 


0.184 


0.024 


0.004 


0.405 


0.135 


0.044 



Table 2: Parameter estimates for generalization I using different values ofv, 5 = 0.5, and A = 0.5 for total 
jump sizes m = 100, 1000, 10000. 









Bias 






RMSE 




V 


Est 


m= 100 


1000 


10000 


m = 100 


1000 


10000 




V 


0.217 


0.079 


-0.022 


0.280 


0.176 


0.124 


0.5 


f 


-0.038 


-0.016 


0.001 


0.072 


0.034 


0.016 




X 


-0.041 


-0.014 


0.006 


0.090 


0.044 


0.030 




V 


0.136 


-0.021 


0.000 


0.222 


0.153 


0.100 


0.6 


f 


-0.026 


0.018 


0.051 


0.067 


0.032 


0.016 




X 


-0.017 


0.005 


0.001 


0.085 


0.041 


0.023 




V 


0.052 


-0.022 


-0.008 


0.188 


0.149 


0.060 


0.7 


f 


-0.015 


0.001 


0.001 


0.071 


0.034 


0.013 




X 


-0.003 


0.007 


0.002 


0.086 


0.040 


0.014 




V 


-0.013 


-0.030 


-0.003 


0.172 


0.129 


0.034 


0.8 


f 


0.004 


0.006 


0.001 


0.073 


0.035 


0.011 




X 


0.007 


0.008 


0.000 


0.086 


0.035 


0.009 




V 


-0.043 


-0.004 


0.000 


0.131 


0.044 


0.013 


0.95 


f 


0.021 


0.002 


0.000 


0.083 


0.028 


0.008 




X 


0.014 


0.002 


0.000 


0.073 


0.022 


0.007 



Table 3: Parameter estimates for generalization II using different values ofv, A = 0.5, and y = 0.5 for total 
jump sizes m = 100, 1000, 10000. 
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